
#simulation study code for paper: Relating latent class assignments to external variables: 
standard errors for correct inference 

nsim <- 500
source("readParam.R")
source("run.template.R")
library(brew)

#simulation conditions for LC entropy: 
alfabeta <- c(1.386294361, 2.197224577)
samplesize <- (c( 'f500', 'f1000', 'f2000')) 
paramsample <- expand.grid(alfabeta=alfabeta, samplesize=samplesize, stringsAsFactors=FALSE)

# options for SE and assignment method: 
SEs <- c("standarderrors", "standarderrors=robust")
assignment <- c('modal', 'prop')


#forloop allconditions in alfabeta 
for(icondition in 1:length(paramsample[,1])) {

    # the brew file for generating data 
    # Here you want to get simulated data, and set stepone analysis to use
    # simulated data set. Loop over nsim
    for (i in 1:nsim){
    envir <- new.env()
    assign("alfabeta", paramsample[icondition, 'alfabeta'], envir=envir)
    assign("samplesize", paramsample[icondition, 'samplesize'], envir=envir)
    run.template("generate.brew", envir=envir, temp.filename.base="generate")

    # brew file for step one analysis: 

        for (iseloop in 1:length(SEs)) {
if (SEs[iseloop] == "standarderrors") {


	envir <- new.env()
	assign("SEs", SEs[iseloop], envir=envir)
		assign("alfabeta", paramsample[icondition, 'alfabeta'], envir=envir) 
	run.template("stepone.brew", envir=envir, temp.filename.base="stepone")



envir <- new.env()
	assign("SEs", SEs[iseloop], envir=envir)
		assign("alfabeta", paramsample[icondition, 'alfabeta'], envir=envir) 
	run.template("onestep.brew", envir=envir, temp.filename.base="onestep")

  
#step 3 - read in D file (containing classifictaion error probabiloities) and create new Dfiles (for the different SE options- that is full
	#file for full correction, file without gradients for 1st order correction, and file with only cl. error logits for no SE correction). 
. 

	# conditions for 3step - inside conditions for step1. 
#loop  over:  - assignment
		# - setype
		# - Dmatrices
	
	    for(istep3 in 1:length(assignment)) {
		
		readParam(paste("Dmatrix", assignment[istep3],
				".txt", sep=""), 1, paste("D", assignment[istep3], ".txt", sep=""))
		readParam(paste("Dmatrix", assignment[istep3], ".txt",
				sep=""), 2, paste("DnoR", assignment[istep3], ".txt", sep=""))

		Dfile <- c(paste("D", assignment[istep3], ".txt", sep=""),
		    paste("DnoR",assignment[istep3], ".txt", sep=""), 
		    paste("Dmatrix",assignment[istep3], ".txt", sep=""))

#for weights
 for (ise2 in 1:length(SEs)) {

		envir <- new.env()
		assign("SEs", SEs[ise2], envir=envir)
			assign("assignment", assignment[istep3], envir=envir)
		

		for (j in 1:3) {

		    assign("Dfile", Dfile[j], envir=envir)
		    run.template("step3.brew", envir=envir, temp.filename.base=
				 paste("c3thirdstep", istep3, j, 
				       assignment[istep3], sep=""))
		 
		}
  #for j loop
	    } #forise2
	}#foristep3 assignment 
} # for if statement SEstep1 one 

# run over condiitons with robust standard errors

if (SEs[iseloop] == "standarderrors=robust") {
#only robust step one and step 3 

	envir <- new.env()
	assign("SEs", SEs[iseloop], envir=envir)
		assign("alfabeta", paramsample[icondition, 'alfabeta'], envir=envir) 
	run.template("stepone.brew", envir=envir, temp.filename.base="stepone")



envir <- new.env()
	assign("SEs", SEs[iseloop], envir=envir)
		assign("alfabeta", paramsample[icondition, 'alfabeta'], envir=envir) 
	run.template("onestep.brew", envir=envir, temp.filename.base="onestep")



#step 3 - read in D file and create new dfiles. 

	# conditions for 3step - inside conditions for step1. 
#loop  over assignment
	#over setype
	# Dmatrices
	
	    for(istep32 in 1:length(assignment)) {
		
		readParam(paste("Dmatrix", assignment[istep32],
				".txt", sep=""), 1, paste("D", assignment[istep32], ".txt", sep=""))
		readParam(paste("Dmatrix", assignment[istep32], ".txt",
				sep=""), 2, paste("DnoR", assignment[istep32], ".txt", sep=""))

		Dfile <- c(paste("D", assignment[istep32], ".txt", sep=""),
		    paste("DnoR",assignment[istep32], ".txt", sep=""), 
		    paste("Dmatrix",assignment[istep32], ".txt", sep=""))

#for weights

		envir <- new.env()
		assign("SEs", SEs[iseloop], envir=envir)
			assign("assignment", assignment[istep32], envir=envir)
		
		for (z in 1:3) {

		    assign("Dfile", Dfile[z], envir=envir)
		    run.template("step3.brew", envir=envir, temp.filename.base=
				 paste("c3thirdstep", istep32, z, 
				       assignment[istep32], sep=""))
		 
		}  #for j loop
	    	}#foristep3 assignment 
} # for if statement SEstep1 one 
}#nsim
}#icondioitn
} # close all end 


#Summarizing  results 

# get a summ over all mean of se, sd, par read in each file separatly, select relevant rows and summarize
# s1: create a file separatly for all conditions with only param and SE v all in one. 
# loop over all files and save to one file: add column for condition and case nr. 

Ds <- c('D', 'DnoR', 'Dmatrix')

library(foreign)

for (i in 1:6) #conditions
	{
	for (j in 1:2)# ses 
		{
 if (j==1) {

		for(k in 1:2)#ses2
			{
 			for (z in 1:length(assignment))
{
				for (y in 1: length(Ds))
					{

    filename2 <- paste ("results",i,j,k,Ds[y],assignment[z],".txt.csv", sep="")
 proba<- read.csv(filename2, header=FALSE, sep="," , 
                   quote="\"",dec=".")

Param <- proba[proba$V5=="Parameters",]
SEs <- proba[proba$V5=="Parameters-SE",]
results <- cbind( 1:length(Param[,1]),i,j,k,Ds[y],assignment[z],Param[,6:13], SEs[,6:13])
write.table(results,"resultsall.dat", sep=" ", append=TRUE, col.names=FALSE)
  }#k
}#if
}#z
}#y

if(j==2) {

 			for (z in 1:length(assignment))
{
				for (y in 1: length(Ds))
					{

    filename2 <- paste ("results",i,j,"2",Ds[y],assignment[z],".txt.csv", sep="")
 proba<- read.csv(filename2, header=FALSE, sep="," , 
                   quote="\"",dec=".")

Param <- proba[proba$V5=="Parameters",]
SEs <- proba[proba$V5=="Parameters-SE",]
results <- cbind( 1:length(Param[,1]),i,j,"2",Ds[y],assignment[z],Param[,6:13], SEs[,6:13])
write.table(results,"resultsall.dat", sep=" ", append=TRUE, col.names=FALSE)
  }#z
}#y
}#if

} #se

}#cond

  

#************************SummarizeONESTEP


#create longfilewith all results 
conFIML <- expand.grid(1:6, 1:2)

for (i in 1:6) #conditions
	{
	for (j in 1:2)# ses 
		{
 
    filename2 <- paste ("resultsFIML",i,j,".csv", sep="")
 proba<- read.csv(filename2, header=FALSE, sep="," , 
                   quote="\"",dec=".")

Param <- proba[proba$V5=="Parameters",]
SEs <- proba[proba$V5=="Parameters-SE",]
results <- cbind( 1:length(Param[,1]),i,j,Param[,6:13], SEs[,6:13])
write.table(results,"resultsFIML.dat", sep=" ", append=TRUE, col.names=FALSE)
  }#k
}#if

#summaryover all conditions for all parameters SD and par estimate 

datFL <- read.table("resultsFIML.dat")

#means separatly per condition 
truevalues <- c( 1.4824, -3.69929, -2, 1,1,0,0,0)

#ge
for (i in 1:6) #conditions
	{
	for (j in 1:2)# ses 
		{
data <- datFL[c(datFL[,3] == i & datFL[,4] == j),]

numdata <- matrix(apply(data,1,function(x){return(as.numeric(x[5:20]))})
,nrow(data),length(5:20),byrow=TRUE)

meanspar <- as.matrix(apply(numdata[,1:8], 2, mean))
meansse <-  as.matrix(apply(numdata[,9:16],2, sd))
sds <- as.matrix(apply(numdata[,1:8],2, sd))

dat <- cbind(i,j, meanspar, sds, meansse)

write.table(dat,"FIMLsumall.dat", sep=" ", append=TRUE, col.names=FALSE)
  }#k
}#if

# for all conditions
	for (j in 1:2)# ses 
		{
data <- datFL[ datFL[,4] == j, ]

numdata <- matrix(apply(data,1,function(x){return(as.numeric(x[5:20]))})
,nrow(data),length(5:20),byrow=TRUE)

meanspar <- as.matrix(apply(numdata[,1:8], 2, mean))
meansse <-  as.matrix(apply(numdata[,9:16],2, sd))
sds <- as.matrix(apply(numdata[,1:8],2, sd))

dat <- cbind(j, meanspar, sds, meansse)

write.table(dat,"FIMLsumTOT.dat", sep=" ", append=TRUE, col.names=FALSE)
  }#k

#get CI as well for step one per condition 
for (i in 1:6) #conditions
	{
	for (j in 1:2)# ses 
		{
data <- datFL[c(datFL[,3] == i & datFL[,4] == j),]

numdata <- matrix(apply(data,1,function(x){return(as.numeric(x[5:20]))})
,nrow(data),length(5:20),byrow=TRUE)

meanspar <- as.matrix(apply(numdata[,1:8], 2, mean))
meansse <-  as.matrix(apply(numdata[,9:16],2, sd))
sds <- as.matrix(apply(numdata[,1:8],2, sd))
Cilower <- as.matrix(numdata[,1:8]-numdata[,9:16])
Ciupper <- as.matrix(numdata[,1:8]+ numdata[,9:16])
#loop over length truevalues
CI <- matrix(NA,499,8)

for(k in 1:499)
{
CI[k,] <- ((Cilower[k,]<=truevalues)+0)*((Ciupper[k,]>=truevalues)+0)
CIrate <- cbind(i,j,t(as.matrix(apply(CI,2,mean))))

}
dat <- cbind(i,j, meanspar, sds, meansse)

write.table(dat,"FIMLsumTOTcompare.dat", sep=" ", append=TRUE, col.names=FALSE)
write.table(CIrate, "covrateFIML.dat", sep=" ", append=TRUE, col.names=FALSE)

}} 





